################################################
################## caculte BLB
r=100
iter.num=500
basis_sd=beta.sd
basis_ci=c0
#####
sd_lst<-list()
ci_lst<-list()
con06<-converge(beta_06,r,basis_sd,basis_ci,not_zero)
con07<-converge(beta_07,r,basis_sd,basis_ci,not_zero)
con08<-converge(beta_08,r,basis_sd,basis_ci,not_zero)
con09<-converge(beta_09,r,basis_sd,basis_ci,not_zero)
conboot<-converge_boot(beta_boot,iter.num,basis_sd,basis_ci,not_zero)
sd_lst[[1]]<-con06$res_sd
sd_lst[[2]]<-con07$res_sd
sd_lst[[3]]<-con08$res_sd
sd_lst[[4]]<-con09$res_sd
sd_lst[[5]]<-conboot$res_sd[-c(1,length(conboot$res_sd))]
ci_lst[[1]]<-con06$res_ci
ci_lst[[2]]<-con07$res_ci
ci_lst[[3]]<-con08$res_ci
ci_lst[[4]]<-con09$res_ci
ci_lst[[5]]<-conboot$res_ci[-c(1,length(conboot$res_ci))]
time_lst<-list()
time_lst[[1]]<-as.vector(time_06$x)
time_lst[[2]]<-as.vector(time_07$x)
time_lst[[3]]<-as.vector(time_08$x)
time_lst[[4]]<-as.vector(time_09$x)
time_lst[[5]]<-as.vector(time_boot[c(2:(iter.num-1)),1])
############
x=sd_lst
y=time_lst
pngname="/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/SD_la_lm.pdf"
my_ylab="SD Relative Deviation"
plotfun(time_lst,sd_lst,pngname,my_ylab=my_ylab,
my_xlim=c(0,400))
######################################
x=ci_lst
y=time_lst
pngname="/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/CI_la_lm.pdf"
my_ylab="CI Relative Deviation"
plotfun(time_lst,ci_lst,pngname,my_ylab=my_ylab,
my_xlim=c(0,400))
################################################
prop_lst=list()
prop_lst[[1]]=apply(beta_06,2,function(x){sum(x!=0)})/nrow(beta_06)
prop_lst[[2]]=apply(beta_07,2,function(x){sum(x!=0)})/nrow(beta_07)
prop_lst[[3]]=apply(beta_08,2,function(x){sum(x!=0)})/nrow(beta_08)
prop_lst[[4]]=apply(beta_09,2,function(x){sum(x!=0)})/nrow(beta_09)
prop_lst[[5]]=apply(beta_boot,2,function(x){sum(x!=0)})/nrow(beta_boot)
pngname="/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/Selection_la_lm.pdf"
seleplot(prop_lst,pngname)
getwd()
setwd(getwd())
seleplot<-function(x,filename){
cols <- gray(0:4/5)
pdf(file=filename, width = 25,height=12, bg = "white")
p=length(x[[1]])
par(mar=c(3,5,2,1),mfrow=c(2,3),cex.axis=2,cex.lab=2,
cex.main=2)
mycol=c(rep("gray",12),rep("white",3),
rep("gray",6),rep("white",4),
rep("gray",8),rep("white",2))
df.bar<-barplot(x[[1]],main="BLBVS-0.6",
ylab="Selection Proportion",col=mycol)
lines(x=df.bar,y=rep(0.5,p))
df.bar<-barplot(x[[2]],main="BLBVS-0.7",
ylab="Selection Proportion",col=mycol)
lines(x=df.bar,y=rep(0.5,p))
df.bar<-barplot(x[[3]],main="BLBVS-0.8",
ylab="Selection Proportion",col=mycol)
lines(x=df.bar,y=rep(0.5,p))
df.bar<-barplot(x[[4]],main="BLBVS-0.9",
ylab="Selection Proportion",col=mycol)
lines(x=df.bar,y=rep(0.5,p))
df.bar<-barplot(x[[5]],main="BootVS",
ylab="Selection Proportion",col=mycol)
lines(x=df.bar,y=rep(0.5,p))
dev.off()
}
plotfun<-function(x,y,filename,my_ylab,my_xlim){
pdf(file=filename, width = 8,height = 6, bg = "white")
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,
cex.lab=2.2,lty=1,lwd=2)
plot(x=x[[1]],y=y[[1]],type="l",
xlab="Time(sec)",ylab=my_ylab,xlim=my_xlim,ylim=c(0,max(y[[5]])*0.5),col="red")
lines(x=x[[2]],y=y[[2]],type="l",col="blue")
lines(x=x[[3]],y=y[[3]],type="l",col="green")
lines(x=x[[4]],y=y[[4]],type="l",col="orange")
lines(x=x[[5]],y=y[[5]],type="l",col="black")
legend(max(my_xlim)-100,max(y[[5]])*0.5,
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8",
"BLBVS-0.9","BootVS"),
lwd=2,pt.cex=1.5,cex=1,
col=c("red","blue","green","orange","black"))
dev.off()
}
converge<-function(beta_blb,r,basis_sd,basis_ci,not_zero){
p=ncol(beta_blb)
s=nrow(beta_blb)/r
std.bag=matrix(0,nrow=s,ncol=p)
CI.bag=matrix(0,nrow=s,ncol=p)
for(i in 1:s){
std.bag[i,]=apply(beta_blb[c(1:r)+(i-1)*r,],2,sd)
CI.bag[i,]=apply(beta_blb[c(1:r)+(i-1)*r,],2,function(x){quantile(x,0.975)-quantile(x,0.025)})
}
######累计平均标准差
std_cum=apply(std.bag,2,cumsum)
std_avg=std_cum
for(i in 1:nrow(std_cum)){
std_avg[i,]=std_cum[i,]/rep(i,ncol(std_cum))
}
###### 累计置信区间宽度
CI_cum=apply(CI.bag,2,cumsum)
CI_avg=CI_cum
for(i in 1:nrow(std_cum)){
CI_avg[i,]=CI_cum[i,]/rep(i,ncol(CI_cum))
}
res_sd=rep(0,nrow(std_avg))
for(i in 1:nrow(std_avg)){
res_sd[i]=mean((abs(std_avg[i,not_zero]-basis_sd[not_zero]))/basis_sd[not_zero])
}
res_ci=rep(0,nrow(CI_avg))
for(i in 1:nrow(CI_avg)){
res_ci[i]=mean((abs(CI_avg[i,not_zero]-basis_ci[not_zero]))/basis_ci[not_zero])
}
return(list("res_sd"=res_sd,"res_ci"=res_ci))
}
converge_boot<-function(beta_boot,r,basis_sd,basis_ci,not_zero){
p=ncol(beta_boot)
std_avg=matrix(0,r,p)
for(i in 2:nrow(beta_boot)){
for(j in 1:ncol(beta_boot)){
std_avg[i,j]=sd(beta_boot[1:i,j])
}
}
CI_avg=matrix(0,r,p)
for(i in 2:nrow(beta_boot)){
for(j in 1:ncol(beta_boot)){
CI_avg[i,j]=quantile(beta_boot[1:i,j],0.975)-quantile(beta_boot[1:i,j],0.025)
}
}
###
res_sd=rep(0,nrow(std_avg))
for(i in 1:nrow(std_avg)){
res_sd[i]=mean((abs(std_avg[i,not_zero]-basis_sd[not_zero]))/basis_sd[not_zero])
}
res_ci=rep(0,nrow(CI_avg))
for(i in 1:nrow(CI_avg)){
res_ci[i]=mean((abs(CI_avg[i,not_zero]-basis_ci[not_zero]))/basis_ci[not_zero])
}
return(list("res_sd"=res_sd,"res_ci"=res_ci))
}
######################################### read data
coef.simu<-read.csv("./beta_hat.csv")
beta_06<-read.csv("./beta_06.csv")
beta_07<-read.csv("./beta_07.csv")
beta_08<-read.csv("./beta_08.csv")
beta_09<-read.csv("./beta_09.csv")
beta_boot<-read.csv("./beta_boot.csv")
colnames(beta_boot)=colnames(beta_09)
time_06<-read.csv("./time_06.csv")
time_07<-read.csv("./time_07.csv")
time_08<-read.csv("./time_08.csv")
time_09<-read.csv("./time_09.csv")
time_boot<-read.csv("./time_boot.csv")
################################################
################## caculate basis
beta.sd=apply(coef.simu,2,sd)
c0=apply(coef.simu,2,function(x){quantile(x,0.975)-quantile(x,0.025)})
not_zero=apply(coef.simu,2,function(x){sum(x!=0)})/nrow(coef.simu)>0.8
################################################
################## caculte BLB
r=100
iter.num=500
basis_sd=beta.sd
basis_ci=c0
#####
sd_lst<-list()
ci_lst<-list()
con06<-converge(beta_06,r,basis_sd,basis_ci,not_zero)
con07<-converge(beta_07,r,basis_sd,basis_ci,not_zero)
con08<-converge(beta_08,r,basis_sd,basis_ci,not_zero)
con09<-converge(beta_09,r,basis_sd,basis_ci,not_zero)
conboot<-converge_boot(beta_boot,iter.num,basis_sd,basis_ci,not_zero)
sd_lst[[1]]<-con06$res_sd
sd_lst[[2]]<-con07$res_sd
sd_lst[[3]]<-con08$res_sd
sd_lst[[4]]<-con09$res_sd
sd_lst[[5]]<-conboot$res_sd[-c(1,length(conboot$res_sd))]
ci_lst[[1]]<-con06$res_ci
ci_lst[[2]]<-con07$res_ci
ci_lst[[3]]<-con08$res_ci
ci_lst[[4]]<-con09$res_ci
ci_lst[[5]]<-conboot$res_ci[-c(1,length(conboot$res_ci))]
time_lst<-list()
time_lst[[1]]<-as.vector(time_06$x)
time_lst[[2]]<-as.vector(time_07$x)
time_lst[[3]]<-as.vector(time_08$x)
time_lst[[4]]<-as.vector(time_09$x)
time_lst[[5]]<-as.vector(time_boot[c(2:(iter.num-1)),1])
############
x=sd_lst
y=time_lst
pngname="/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/fig/SD_la_lm.pdf"
my_ylab="SD Relative Deviation"
plotfun(time_lst,sd_lst,pngname,my_ylab=my_ylab,
my_xlim=c(0,400))
######################################
x=ci_lst
y=time_lst
pngname="/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/fig/CI_la_lm.pdf"
my_ylab="CI Relative Deviation"
plotfun(time_lst,ci_lst,pngname,my_ylab=my_ylab,
my_xlim=c(0,400))
################################################
prop_lst=list()
prop_lst[[1]]=apply(beta_06,2,function(x){sum(x!=0)})/nrow(beta_06)
prop_lst[[2]]=apply(beta_07,2,function(x){sum(x!=0)})/nrow(beta_07)
prop_lst[[3]]=apply(beta_08,2,function(x){sum(x!=0)})/nrow(beta_08)
prop_lst[[4]]=apply(beta_09,2,function(x){sum(x!=0)})/nrow(beta_09)
prop_lst[[5]]=apply(beta_boot,2,function(x){sum(x!=0)})/nrow(beta_boot)
pngname="/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/fig/Selection_la_lm.pdf"
seleplot(prop_lst,pngname)
x<-prop_lst
filename=pngname
cols <- gray(0:4/5)
pdf(file=filename, width = 25,height=12, bg = "white")
p=length(x[[1]])
par(mar=c(3,5,2,1),mfrow=c(2,3),cex.axis=2,cex.lab=2,
cex.main=2)
mycol=c(rep("gray",12),rep("white",3),
rep("gray",6),rep("white",4),
rep("gray",8),rep("white",2))
df.bar<-barplot(x[[1]],main="BLBVS-0.6",
ylab="Selection Proportion",col=mycol)
lines(x=df.bar,y=rep(0.5,p))
df.bar<-barplot(x[[2]],main="BLBVS-0.7",
ylab="Selection Proportion",col=mycol)
lines(x=df.bar,y=rep(0.5,p))
df.bar<-barplot(x[[3]],main="BLBVS-0.8",
ylab="Selection Proportion",col=mycol)
lines(x=df.bar,y=rep(0.5,p))
df.bar<-barplot(x[[4]],main="BLBVS-0.9",
ylab="Selection Proportion",col=mycol)
lines(x=df.bar,y=rep(0.5,p))
df.bar<-barplot(x[[5]],main="BootVS",
ylab="Selection Proportion",col=mycol)
lines(x=df.bar,y=rep(0.5,p))
dev.off()
############
x=sd_lst
y=time_lst
pngname="/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/fig/SD_la_lm.pdf"
my_ylab="SD Relative Deviation"
plotfun(time_lst,sd_lst,pngname,my_ylab=my_ylab,
my_xlim=c(0,400))
setwd(getwd())
beta_06<-read.csv("./beta_06.csv")
beta_07<-read.csv("./beta_07.csv")
beta_08<-read.csv("./beta_08.csv")
beta_09<-read.csv("./beta_09.csv")
beta_boot<-read.csv("./beta_boot.csv")
p<-ncol(beta_06)
prop_06<-sort(apply(beta_06,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_07<-sort(apply(beta_07,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_08<-sort(apply(beta_08,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_09<-sort(apply(beta_09,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_boot<-sort(apply(beta_boot,2,function(x)sum(x!=0)/length(x)),decreasing=T)
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,cex.lab=2.2,lty=1,lwd=2,family='STKaiti')
plot(c(1:35),prop_06,
type="o",ylab="比例",xlab=" ",
lty="longdash",ylim=c(0,1),pch=1)
lines(c(1:35),prop_07,
type="o",lty="dashed",pch=2)
lines(c(1:35),prop_08,
type="o",lty="dotted",pch=3)
lines(c(1:35),prop_09,
type="o",lty="dotdash",pch=4)
lines(c(1:35),prop_boot,
type="o",lty="solid",pch=5)
legend(1,0.9,
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8","BLBVS-0.9","BootVS"),
lty=c("longdash","dashed","dotted","dotdash","solid"),
lwd=2,pt.cex=1.4,cex=1.7,pch=c(1:5))
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,cex.lab=2.2,lty=1,lwd=2)
plot(c(1:35),prop_06,
type="o",ylab="Selection Proportion",xlab=" ",
col="red",ylim=c(0,1),pch=1)
lines(c(1:35),prop_07,
type="o",col="blue",pch=2)
lines(c(1:35),prop_08,
type="o",col="green",pch=3)
lines(c(1:35),prop_09,
type="o",col="orange",pch=4)
lines(c(1:35),prop_boot,
type="o",col="black",pch=5)
legend(1,0.9,
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8","BLBVS-0.9","BootVS"),
col=c("red","blue","green","orange","black"),
lwd=2,pt.cex=1.4,cex=1.7,pch=c(1:5))
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,cex.lab=2.2,lty=1,lwd=2)
plot(c(1:35),prop_06,
type="o",ylab="Selection Proportion",xlab=" ",
col="red",ylim=c(0,1),pch=1)
lines(c(1:35),prop_07,
type="o",col="blue",pch=1)
lines(c(1:35),prop_08,
type="o",col="green",pch=1)
lines(c(1:35),prop_09,
type="o",col="orange",pch=1)
lines(c(1:35),prop_boot,
type="o",col="black",pch=1)
legend(1,0.9,
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8","BLBVS-0.9","BootVS"),
col=c("red","blue","green","orange","black"),
lwd=2,pt.cex=1.4,cex=1.7,pch=c(1:5))
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,cex.lab=2.2,lty=1,lwd=2)
plot(c(1:35),prop_06,
type="o",ylab="Selection Proportion",xlab=" ",
col="red",ylim=c(0,1),pch=1)
lines(c(1:35),prop_07,
type="o",col="blue",pch=1)
lines(c(1:35),prop_08,
type="o",col="green",pch=1)
lines(c(1:35),prop_09,
type="o",col="orange",pch=1)
lines(c(1:35),prop_boot,lty="dashed"
type="o",col="black",pch=1)
legend(1,0.9,
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8","BLBVS-0.9","BootVS"),
col=c("red","blue","green","orange","black"),
lwd=2,pt.cex=1.4,cex=1.7,pch=c(1:5))
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,cex.lab=2.2,lty=1,lwd=2)
plot(c(1:35),prop_06,
type="o",ylab="Selection Proportion",xlab=" ",
col="red",ylim=c(0,1),pch=1)
lines(c(1:35),prop_07,
type="o",col="blue",pch=1)
lines(c(1:35),prop_08,
type="o",col="green",pch=1)
lines(c(1:35),prop_09,
type="o",col="orange",pch=1)
lines(c(1:35),prop_boot,lty="dashed",
type="o",col="black",pch=1)
legend(1,0.9,
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8","BLBVS-0.9","BootVS"),
col=c("red","blue","green","orange","black"),
lwd=2,pt.cex=1.4,cex=1.7,pch=c(1:5))
png(filename="/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/fig/Cutoff_la_lm.PNG", width = 25,height = 18, units ="cm", bg = "white", res = 600)
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,cex.lab=2.2,lty=1,lwd=2)
plot(c(1:35),prop_06,
type="o",ylab="Selection Proportion",xlab=" ",
col="red",ylim=c(0,1),pch=1)
lines(c(1:35),prop_07,
type="o",col="blue",pch=1)
lines(c(1:35),prop_08,
type="o",col="green",pch=1)
lines(c(1:35),prop_09,
type="o",col="orange",pch=1)
lines(c(1:35),prop_boot,
type="o",col="black",pch=1)
legend(1,0.9,
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8","BLBVS-0.9","BootVS"),
col=c("red","blue","green","orange","black"),
lwd=2,pt.cex=1.4,cex=1.7,pch=c(1:5))
dev.off()
pdf(file="/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/fig/Cutoff_la_lm.pdf", width = 25,height = 18, bg = "white")
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,cex.lab=2.2,lty=1,lwd=2)
plot(c(1:35),prop_06,
type="o",ylab="Selection Proportion",xlab=" ",
col="red",ylim=c(0,1),pch=1)
lines(c(1:35),prop_07,
type="o",col="blue",pch=1)
lines(c(1:35),prop_08,
type="o",col="green",pch=1)
lines(c(1:35),prop_09,
type="o",col="orange",pch=1)
lines(c(1:35),prop_boot,
type="o",col="black",pch=1)
legend(1,0.9,
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8","BLBVS-0.9","BootVS"),
col=c("red","blue","green","orange","black"),
lwd=2,pt.cex=1.4,cex=1.7,pch=c(1:5))
dev.off()
setwd(getwd())
beta_06<-read.csv("./beta_06.csv")
beta_07<-read.csv("./beta_07.csv")
beta_08<-read.csv("./beta_08.csv")
beta_09<-read.csv("./beta_09.csv")
beta_boot<-read.csv("./beta_boot.csv")
p<-ncol(beta_06)
prop_06<-sort(apply(beta_06,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_07<-sort(apply(beta_07,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_08<-sort(apply(beta_08,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_09<-sort(apply(beta_09,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_boot<-sort(apply(beta_boot,2,function(x)sum(x!=0)/length(x)),decreasing=T)
####选择cut-off
pdf(file="/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/fig/Cutoff_la_lm.pdf", width =8,height = 6, bg = "white")
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,cex.lab=2.2,lty=1,lwd=2)
plot(c(1:35),prop_06,
type="o",ylab="Selection Proportion",xlab=" ",
col="red",ylim=c(0,1),pch=1)
lines(c(1:35),prop_07,
type="o",col="blue",pch=1)
lines(c(1:35),prop_08,
type="o",col="green",pch=1)
lines(c(1:35),prop_09,
type="o",col="orange",pch=1)
lines(c(1:35),prop_boot,
type="o",col="black",pch=1)
legend(1,0.9,
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8","BLBVS-0.9","BootVS"),
col=c("red","blue","green","orange","black"),
lwd=2,pt.cex=1.4,cex=1.7,pch=c(1:5))
dev.off()
pdf(file="/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/fig/Cutoff_la_lm.pdf", width =8,height = 6, bg = "white")
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,cex.lab=2.2,lty=1,lwd=2)
plot(c(1:35),prop_06,
type="o",ylab="Selection Proportion",xlab=" ",
col="red",ylim=c(0,1),pch=1)
lines(c(1:35),prop_07,
type="o",col="blue",pch=1)
lines(c(1:35),prop_08,
type="o",col="green",pch=1)
lines(c(1:35),prop_09,
type="o",col="orange",pch=1)
lines(c(1:35),prop_boot,
type="o",col="black",pch=1)
legend(1,0.9,
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8","BLBVS-0.9","BootVS"),
col=c("red","blue","green","orange","black"),
lwd=2,pt.cex=1.4,cex=1.7)
dev.off()
pdf(file="/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/fig/Cutoff_la_lm.pdf", width =8,height = 6, bg = "white")
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,cex.lab=2.2,lty=1,lwd=2)
plot(c(1:35),prop_06,
type="o",ylab="Selection Proportion",xlab=" ",
col="red",ylim=c(0,1),pch=1)
lines(c(1:35),prop_07,
type="o",col="blue",pch=1)
lines(c(1:35),prop_08,
type="o",col="green",pch=1)
lines(c(1:35),prop_09,
type="o",col="orange",pch=1)
lines(c(1:35),prop_boot,
type="o",col="black",pch=1)
legend(1,0.9,
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8","BLBVS-0.9","BootVS"),
col=c("red","blue","green","orange","black"),
lwd=2,pt.cex=1.4,cex=1.2)
dev.off()
setwd(getwd())
beta_06<-read.csv("./beta_06.csv")
beta_07<-read.csv("./beta_07.csv")
beta_08<-read.csv("./beta_08.csv")
beta_09<-read.csv("./beta_09.csv")
beta_boot<-read.csv("./beta_boot.csv")
p<-ncol(beta_06)
prop_06<-sort(apply(beta_06,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_07<-sort(apply(beta_07,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_08<-sort(apply(beta_08,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_09<-sort(apply(beta_09,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_boot<-sort(apply(beta_boot,2,function(x)sum(x!=0)/length(x)),decreasing=T)
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,cex.lab=2.2,lty=1,lwd=2)
plot(c(1:35),prop_06,
type="o",ylab="Selection Proportion",xlab="Number of Variables",
col="red",ylim=c(0,1),pch=1)
lines(c(1:35),prop_07,
type="o",col="blue",pch=1)
lines(c(1:35),prop_08,
type="o",col="green",pch=1)
lines(c(1:35),prop_09,
type="o",col="orange",pch=1)
lines(c(1:35),prop_boot,
type="o",col="black",pch=1)
legend(1,0.9,
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8","BLBVS-0.9","BootVS"),
col=c("red","blue","green","orange","black"),
lwd=2,pt.cex=1.4,cex=1.2)
setwd(getwd())
beta_06<-read.csv("./beta_06.csv")
beta_07<-read.csv("./beta_07.csv")
beta_08<-read.csv("./beta_08.csv")
beta_09<-read.csv("./beta_09.csv")
beta_boot<-read.csv("./beta_boot.csv")
p<-ncol(beta_06)
prop_06<-sort(apply(beta_06,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_07<-sort(apply(beta_07,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_08<-sort(apply(beta_08,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_09<-sort(apply(beta_09,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_boot<-sort(apply(beta_boot,2,function(x)sum(x!=0)/length(x)),decreasing=T)
####选择cut-off
pdf(file="/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/fig/Cutoff_la_lm.pdf", width =8,height = 6, bg = "white")
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,cex.lab=2.2,lty=1,lwd=2)
plot(c(1:35),prop_06,
type="o",ylab="Selection Proportion",xlab="Number of Variables",
col="red",ylim=c(0,1),pch=1)
lines(c(1:35),prop_07,
type="o",col="blue",pch=1)
lines(c(1:35),prop_08,
type="o",col="green",pch=1)
lines(c(1:35),prop_09,
type="o",col="orange",pch=1)
lines(c(1:35),prop_boot,
type="o",col="black",pch=1)
legend(1,0.9,
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8","BLBVS-0.9","BootVS"),
col=c("red","blue","green","orange","black"),
lwd=2,pt.cex=1.4,cex=1.2)
dev.off()
